clc; clear; close all;
% Lyapunov 指数随 memristor 初始条件 z0 的变化

% 参数设置
mu = 3.5;
x0_init = 0.1; 
y0_init = 0;
z0_init = -1;        % 固定 (x0,y0,z0) = (0.1, 0, -1)

nk = 300;            % k 的采样点数
k_vals = linspace(0, 0.2, nk);  % 扫描 k 的区间

N_transient = 500;   % 瞬态迭代步数
N_LEs = 500;         % 用于计算 Lyapunov 指数的迭代步数

% 初始化存储 Lyapunov 指数的矩阵
LEs_array = zeros(3, nk);

% 扫描 k 并计算 Lyapunov 指数
for i = 1:nk
    k = k_vals(i);
    
    % 1) 瞬态迭代（去除初始影响）
    state = [x0_init, y0_init, z0_init];
    for n = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    
    % 2) 从瞬态后的状态开始计算 Lyapunov 指数
    LEs_result = LEs(state, mu, k, N_LEs);  % 这里 LEs 返回一个 3x1 数组
    LEs_array(:, i) = LEs_result(:);
end

figure;
hold on;
plot(k_vals, LEs_array(1,:), 'r-', 'LineWidth',1);  % LE1
plot(k_vals, LEs_array(2,:), 'b-', 'LineWidth',1);  % LE2
plot(k_vals, LEs_array(3,:), 'm-', 'LineWidth',1);  % LE3
xlabel('k');
ylabel('Lyapunov Exponents');
title('Lyapunov 指数谱与 k 的关系');
legend('LE1','LE2','LE3','Location','best');
grid on;
